% EXPLAINER: 

%% Generates:
% Table 4
% Tables A7-A10

%% File is structured as follows: 
% 0) import data inc set specification for preferences
% 1) bootstrap preference differences
% 2) Bootstrap belief differences

%% 0) import data inc set specification for preferences
clear; clc;
close all;
rng(88);

pref_spec=2; % set preference specification (just use main one, 2)

% load main estimates
%----------------------
load('../Output/main_estimates.mat')

% import data - ex ante
%-----------------------------
[xadata] = import_and_setup_exante;

% import and setup expost data
%-----------------------------
[xpdata,n_prefparms,n_groomchars] = import_and_setup_expost(pref_spec);

% set up states
% columns are: age, completed education, in school, current education (if
% in school)
%-----------------------------
setup_states;


%% 1) bootstrap preference differences

    N=500;
    p_starting_vals=estimates.pref;
     lb=estimates.pref-5;
     ub=estimates.pref+5;

    % fix parameters that we are not testing heterogeneity in
    ub(1:8)=estimates.pref(1:8)
    lb(1:8)=estimates.pref(1:8)
    p_starting_vals(1:8)=estimates.pref(1:8)
    ub(end)=estimates.pref(end)
    lb(end)=estimates.pref(end)
   
    options_f = optimoptions('fmincon', 'UseParallel', false, 'MaxIterations', 1e6, 'CheckGradients', false, 'MaxFunEvals', 1000000, ...
    'SpecifyObjectiveGradient', false, 'OptimalityTolerance', 1e-6, 'SpecifyConstraintGradient',false, 'FiniteDifferenceType',  'forward');

    hetnames={'inschool16', 'married17', 'mothers_ed', 'fathers_ed', 'scst', 'dirtfloor', 'married_child', 'older', 'treat1', 'treat2'}
    for i=1:length(hetnames)
        hetname=hetnames{i}
        
            hello= sprintf('../Data/created_data/xp_w_%s0.csv',hetname);
            hw0 = csvread(hello, 0,0);
            hello= sprintf('../Data/created_data/xp_w_%s1.csv',hetname);
            hw1 = csvread(hello, 0,0); 
    
            % bs
            bs1=ones(size(p_starting_vals',1),N);
            bs0=ones(size(p_starting_vals',1),N);

            parfor n=0:N 
                hetname
                n
                w=hw1(:,n+1);
                [temp,ll,~,~,lambda,grad,H] = fmincon(@(alpha0) pref_mom(xpdata,n_prefparms,alpha0,w), ...
                            estimates.pref, [], [], [], [],lb,ub, [], options_f);
                bs1(:,n+1)=[temp'];
            end

            parfor n=0:N 
                hetname
                n
                w=hw0(:,n+1);
                [temp,ll,~,~,lambda,grad,H] = fmincon(@(alpha0) pref_mom(xpdata,n_prefparms,alpha0,w), ...
                            estimates.pref, [], [], [], [],lb,ub, [], options_f);
                bs0(:,n+1)=[temp'];
            end
                        
           se0= std(bs0(:,2:end),0,2);
           se1= std(bs1(:,2:end),0,2);
            
           real0=bs0(:,1);
           real1=bs1(:,1);
           
            % Table of pref differences
            diff=real1-real0;
            diff_se=(se0.^2+se1.^2).^0.5
            t=diff./diff_se;
            pref_results_het=[real1 se1 real0 se0 diff diff_se t]
            colNames = {'coef1','se1','coef0','se0','difference', 'diff-se','t'}
            sTable = array2table(pref_results_het,'VariableNames',colNames) %,'RowNames',rowNames) 
            
            hello= sprintf('../Output_New/hettable_%s.csv',hetname);
            hetname
            writetable(sTable,hello)
            
            hello= sprintf('../Output_New/bs_results/bs_%s0.csv',hetname);
            csvwrite(hello,bs0);
            
            hello= sprintf('../Output_New/bs_results/bs_%s0',hetname);
            save(hello,'bs0')

            hello= sprintf('../Output_New/bs_results/bs_%s1.csv',hetname);
            csvwrite(hello,bs1);
            
            hello= sprintf('../Output_New/bs_results/bs_%s1',hetname);
            save(hello,'bs1')

    end

%% 2) Bootstrap belief differences
load('../Output/main_estimates.mat')

clear fixedparms


matvalues= ... % starting, lb, ub
            [   estimates.beliefs(1),    estimates.beliefs(1),      estimates.beliefs(1);           ... % sigma - constant
                0,      0,      0;      ... % sigma - het in age
                -3,     -8,      1;     ... % beliefs - constant
                0,      -1,      1;     ... % beliefs - age
                1.2,    -2,      3;     ... % beliefs - ed
                0,    0,      0;        ... % beliefs - ed^2
                2,      -1,   5;        ... % beliefs - college
                -0.1,      -0.5,   0.5; ... % beliefs - age x ed
                estimates.beliefs(end-4),     estimates.beliefs(end-4),     estimates.beliefs(end-4); ...
                estimates.beliefs(end-3),     estimates.beliefs(end-3),     estimates.beliefs(end-3); ...
                estimates.beliefs(end-2),     estimates.beliefs(end-2),     estimates.beliefs(end-2); ...
                estimates.beliefs(end-1),     estimates.beliefs(end-1),     estimates.beliefs(end-1); ...  % weight on u0 in VTp1
                estimates.beliefs(end),   estimates.beliefs(end),   estimates.beliefs(end)]                % inattentive share


        starting_vals=matvalues(:,1)';
        lb=matvalues(:,2)';
        ub=matvalues(:,3)';

        n_beliefparms=size(lb,2)-6;
        nvars=size(lb,2);
        rng(1);
        
        N=50
   
        fixedparms.n_prefparms=n_prefparms
        fixedparms.n_beliefparms=n_beliefparms
        fixedparms.n_groomchars=n_groomchars
 
        fixedparms.disc=0.95;
        
        spec=2;

        
    hetnames={'inschool16', 'married17', 'mothers_ed', 'fathers_ed', 'scst', 'dirtfloor', 'married_child', 'older', 'treat1', 'treat2'}
    for i=1:length(hetnames)
        hetname=hetnames{i}

        % open weights
            hello= sprintf('../Data/created_data/xa_w_%s0.csv',hetname);
            hw0 = csvread(hello, 0,0);
            hello= sprintf('../Data/created_data/xa_w_%s1.csv',hetname);
            hw1 = csvread(hello, 0,0);  

            % open preference estimates 
            temp=load('../Output/bs_results/bs_prefs.mat');
            bs_prefs=temp.bs_prefs;
            hello= sprintf('../Output/bs_results/bs_%s0.mat',hetname);
            temp=load(hello);
            bs_prefs0=temp.bs0;
            hello= sprintf('../Output/bs_results/bs_%s1.mat',hetname);
            temp=load(hello);
            bs_prefs1=temp.bs1;

            % bs
            bs1=ones(size(starting_vals',1),N);
            bs0=ones(size(starting_vals',1),N);
    
               
            xadata0.G_exante=xadata.G_exante(hw0(:,1)>0,:);
            xadata0.age=xadata.age(hw0(:,1)>0,:);
            xadata0.ed=xadata.ed(hw0(:,1)>0,:);
            xadata0.stillinschool=xadata.stillinschool(hw0(:,1)>0,:);
            xadata0.chosen_marry=xadata.chosen_marry(hw0(:,1)>0,:);
            xadata0.chosen_school=xadata.chosen_school(hw0(:,1)>0,:);
            xadata0.good_girl=xadata.good_girl(hw0(:,1)>0,:);
            xadata0.like_school=xadata.like_school(hw0(:,1)>0,:);
            hw0=hw0(hw0(:,1)>0,:);
            
            xadata1.G_exante=xadata.G_exante(hw1(:,1)>0,:);
            xadata1.age=xadata.age(hw1(:,1)>0,:);
            xadata1.ed=xadata.ed(hw1(:,1)>0,:);
            xadata1.stillinschool=xadata.stillinschool(hw1(:,1)>0,:);
            xadata1.chosen_marry=xadata.chosen_marry(hw1(:,1)>0,:);
            xadata1.chosen_school=xadata.chosen_school(hw1(:,1)>0,:);
            xadata1.good_girl=xadata.good_girl(hw1(:,1)>0,:);
            xadata1.like_school=xadata.like_school(hw1(:,1)>0,:);
            hw1=hw1(hw1(:,1)>0,:);
            
            options = optimoptions('particleswarm','InitialSwarmMatrix',starting_vals,'UseParallel',true,'FunctionTolerance',1e-4,'OutputFcn',{@pigraph},'SwarmSize',100,'HybridFcn',...
                            @fmincon,'Display','iter')

          
            % 1 
            for n=0:N
                n
                w=hw1(:,n+1);
                fixedparms.prefparms=bs_prefs1(:,n+1)';
                fun_bs=@(alpha0) aggregate_beliefs(states,xadata1,fixedparms,spec,w,alpha0); 
                rng(1);
               [temp,f] = particleswarm(fun_bs,nvars,lb,ub,options)
                bs1(:,n+1)=temp';
            end 
               

    
            % 0 
            for n=0:N
                n
                w=hw0(:,n+1);
                fixedparms.prefparms=bs_prefs0(:,n+1)';
                fun_bs=@(alpha0) aggregate_beliefs(states,xadata0,fixedparms,spec,w,alpha0); 
                rng(1);
                [temp,f] = particleswarm(fun_bs,nvars,lb,ub,options)
                bs0(:,n+1)=temp';
            end 
            
           se0= std(bs0(:,2:end),0,2);
           se1= std(bs1(:,2:end),0,2);
            
           real0=bs0(:,1);
           real1=bs1(:,1);
           

            % Table of belief differences
            diff=real1-real0;
            diff_se=(se0.^2+se1.^2).^0.5
            t=diff./diff_se;
            belief_results_het=[real1 se1 real0 se0 diff diff_se t]
            colNames = {'coef1','se1','coef0','se0','difference', 'diff-se','t'}
            sTable = array2table(belief_results_het,'VariableNames',colNames) 
            
            
            hello= sprintf('../Output_New/hettable_b_%s.csv',hetname);
            hetname
            writetable(sTable,hello)
            
            hello= sprintf('../Output_New/bs_results/bs_b_%s0.csv',hetname);
            csvwrite(hello,bs0);
            
            hello= sprintf('../Output_New/bs_results/bs_b_%s0',hetname);
            save(hello,'bs0')

            hello= sprintf('../Output_New/bs_results/bs_b_%s1.csv',hetname);
            csvwrite(hello,bs1);
            
            hello= sprintf('../Output_New/bs_results/bs_b_%s1',hetname);
            save(hello,'bs1')
            
   end
  